Learning Objectives

  • Introduce map outlines available in the maps package.
  • Learn to set and transform coordinate systems using the sf package.
  • Reshape and plot spatial data using the tidyverse.
  • Plot raster data using the terrainr and raster packages.
  • Add scale bars and north arrows using the ggspatial package.
  • Next week: Apply skills to create a map of field sites or sampling locations.

Introduction

Today we will cover the basics of creating maps in R. By the end of this workshop, you will be able to create a map and add labeled points, raster data, and shape files. Our introduction to mapping will, in many ways, be an extension of data manipulation and plotting within the tidyverse.

Getting started

Open your Environmental Data Analysis project, and create a new script named “A10_Mapping”. Download the mapfiles folder from Brightspace, and extract the files into a new folder named “mapfiles” in your project folder. Remember to back up your project folder after completing the workshop today.

Packages

Load the packages you will need today. Computers 1-20 in Hudson 032 should have all of these packages installed already. If you are working on your personal computers, you will need to install the new mapping packages.

#Load packages
library(tidyverse)
library(maps)
library(sf)
library(terrainr)
library(raster)
library(ggspatial)

Retrieving basic map data

The first thing we will do is download a map of the continental United States. The maps package contains a file named “state” that contains georeferenced polygons for all of the lower 48 states. We will use the map() function to import these data and assign them to the object state. The fill = T argument specifies that you would like each state to be stored as a separate polygon rather than a series of lines.

# Get states
state <- map("state", fill = T)

This function should return an outline of the continental United States. Note: The package purrr (part of the tidyverse) and the map package both contain the function map(). This may cause an error when you attempt to retrieve the map data. When you are working with packages that contain functions of the same name, you can specify which package you mean by adding the package name ahead of the function using the syntax package::function(). So for the example above, you would specify that you want to use the map function from the maps package using maps::map("state", fill = T).

Transforming data into a spatial object

If you look at your state object in the environment window, you will see that the function has created a list with four vectors. These vectors contain x and y coordinates (longitude and latitude), the bounds of the area, and the names of the states. We would like to turn this list into a spatial object that we can use to visualize and analyze relationships between locations. The sf package allows us to do just that. We will use the st_as_sf() to transform our state map object into a spatial data frame.

# Create sf object
statesf <- st_as_sf(state)

Now if you take a look at the statesf object in the environment window, you will see that the function has created a data frame with each state as a separate row in the ID column (including the District of Columbia) and a geom column that contains a list specifying the x-y coordinates of the outline of the state polygon. Run your object name to see more information that is stored in this object.

statesf
## Simple feature collection with 49 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.6813 ymin: 25.12993 xmax: -67.00742 ymax: 49.38323
## Geodetic CRS:  +proj=longlat +ellps=clrk66 +no_defs +type=crs
## First 10 features:
##                                        ID                           geom
## alabama                           alabama MULTIPOLYGON (((-87.46201 3...
## arizona                           arizona MULTIPOLYGON (((-114.6374 3...
## arkansas                         arkansas MULTIPOLYGON (((-94.05103 3...
## california                     california MULTIPOLYGON (((-120.006 42...
## colorado                         colorado MULTIPOLYGON (((-102.0552 4...
## connecticut                   connecticut MULTIPOLYGON (((-73.49902 4...
## delaware                         delaware MULTIPOLYGON (((-75.80231 3...
## district of columbia district of columbia MULTIPOLYGON (((-77.13731 3...
## florida                           florida MULTIPOLYGON (((-85.01548 3...
## georgia                           georgia MULTIPOLYGON (((-80.89018 3...

Note that the range from the maps object appears as the “Bounding box” in the sf object, the geometry appears as a MULTIPOLYGON, and the coordinate system of the map appears as the Geodetic CRS. In this case, the coordinate reference system (CRS) for our map is WGS 84, which is commonly used in North America.

Once you have created an sf object, you can use tidyverse functions to reshape and plot the data. Let’s use the ggplot() function to create a quick map. We will use ggplot() to create the plot object and identify our data frame, and then add a geom_sf() layer to add our sf object geometry to the map.

ggplot(statesf) +
  geom_sf()

Note You may notice that your degree symbols are not plotting correctly. If so, this is a graphics issue. You can fix this issue in RStudio by selecting Tools -> Global Options. This will open a window. Choose the Graphics tab at the top of the window, and change the Backend to “Cairo.”

Change plot aesthetics

All of the arguments in ggplot can be added to customize your plot. You can change the color of the line and the fill of the polygons, or the width of the outline. You can change the theme using the theme_minimal() or the theme_bw() functions we have already learned. A common theme for mapping is theme_void, which removes all elements except for the map layers.

ggplot(statesf) +
  geom_sf(color = "#283044", fill = "#74D3AE", size = 0.5) +
  theme_void()

Note that fixed modifiers are included outside of a aes() mapping function. If you wished to add a modifier based on a variable, that would need to be included inside the aes() function. For example, if I wanted to plot every state as a different color, I could add the argument aes(fill = ID) to the geom_sf() layer.

ggplot(statesf) +
  geom_sf(aes(fill = ID)) +
  theme_void()

Adding point data

You can also easily plot points on the map based on coordinates in a tibble or data frame. When you add a new layer, in this case a geom_point() and a geom_text() layer, you need only specify the data frame where the information is located using data = and then specify your aesthetics. Let’s add points to show the positions of SUNY Plattsburgh and Washington D.C. on our map.

First we will create a tibble with names and coordinates for the two locations, then we will add points and labels to our map.

# Create table with coordinates and names
labs <- tibble(
  long = c(-73.464279, -77.015541),
  lat = c(44.692614, 38.892753),
  names = c("SUNY Plattsburgh", "Washington D.C."))

# Add point and text layers to map
ggplot(statesf) +
  geom_sf(color = "#F5DBCB", fill = "#74D3AE", size = 0.5) +
  geom_point(data = labs, aes(x = long, y = lat), 
             color = "#283044", size = 5) +
  geom_text(data = labs, aes(x = long, y = lat, label = names),
            hjust = 1, nudge_x = -1, color = "#283044")

Notice that I added the arguments hjust and nudge_x to orient the labels so that they are readable. If you are able to install packages, the package ggrepel has a function geom_label_repel() that you can add in place of geom_text(). This function will automatically space labels so that they are readable.

Challenge 1

Create a plot of the United States with every state plotted as a different color.
Hint: Look at the statesf object to determine which variable to use. Think about where you need to add the fill = argument to change fill based on a variable in the statesf data frame.
Change the color that is used to outline the states.
Consult the ggplot2 cheat sheet or Google to figure out how to exclude the legend from appearing on your map.
When you are satisfied with your map, use ggsave("A10_C1.jpg") to export your map image.
Upload your plot to Challenge 1 on Brightspace.

Answer 
ggplot(statesf) +
  geom_sf(aes(fill = ID), color = "white") +
  theme_minimal() +
  theme(legend.position = "none") 

ggsave("figures/A10_C1.jpg")
## Saving 7 x 5 in image

Plot a subset of states

Because the sf object is stored as a data frame, we can use the tidy tools of dplyr to retain just certain parts of it. For example, we can select just NY and VT and then plot them.

# Create data frame for just NY and VT
ny_vt <- statesf %>% 
  filter(ID %in% c("new york", "vermont"))
  
#Plot only those two states
ggplot(ny_vt) +
  geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
  theme_minimal()

Adding shape files

So far, we’ve added point and text data to our base map by adding layers to our ggplot and specifying a data object. We can use a similar approach to layer additional spatial data onto our map. In the next section, we will import a shape file for Lake Champlain as an sf object, transform it into the same coordinate reference system (CRS) as our current layers, and then plot it using the geom_sf() function.

A shape file is a data storage format that stores the location, shape, and various attributes of spatial features. A shape file is stored as a set of related files (i.e., cpg, dbf, prj, shp, and shx) with the same file name prefix. All of these files are needed to plot the geometry appropriately and should be stored in the same folder in your file directory.

Let’s import our shape file for Lake Champlain from the mapfiles folder using the function st_read() in the sf package. The syntax for this function is st_read("folder name", "file prefix").

# Import shape file data from mapfiles folder
lc <- st_read("mapfiles", "VT_Lake_Champlain")
## Reading layer `VT_Lake_Champlain' from data source 
##   `C:\Users\malld001\Dropbox\PlattsburghCourses\ENV422_EnvDataAnalysis_Spring2024\A10Mapping\mapfiles' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 420142.5 ymin: 115042.3 xmax: 454613.1 ymax: 287432.3
## Projected CRS: NAD83 / Vermont

Take a look at the attributes for the sf object you have just imported. Note that the Projected CRS differs from the “WGS 84” system on the maps we have made up to this point. We should transform the coordinate system of our shape file if we want the layers to match up properly.

Coordinate Reference Systems

When we plot a map in two dimensions, we are necessarily translating 3D information into a 2D picture. This is easier for small areas than large areas, but maps are always a projection of a 3D reality. R can read more than 6,000 different projections and transform them to a different CRS. Yes there are that many! Fortunately, within a region, only ~3-5 are typically used. Most spatial data are stored with this information, or the information appears in the meta-data. For example, our Lake Champlain shape file specifies the CRS as NAD83 / Vermont.

In some cases, you may import (or be given) a file with no CRS associated with it. You will need to figure out the projection system used in order to transform it to the projection you are using. Some of the most commonly used systems in North America are listed in the table below.
Table 10.1 Coordinate Reference Systems (CRS) and associated European Petroleum Survey Group (EPSG) codes commonly used in North America
CRS EPSG Note
WGS84 4326 Commonly used by organizations that provide GIS data for the entire globe or many countries. CFS used by Google Earth.
NAD83 4269 Most commonly used by US federal agencies
NAD27 4267 Old version of NAD83
Mercator 3857 Tiles from Google Maps, Open Street Maps, Stamen Maps
UTM, Zone 10 32610 See map below
UTM, Zone 11 32611 See map below
UTM, Zone 12 32612 See map below
UTM, Zone 13 32613 See map below
UTM, Zone 14 32614 See map below
UTM, Zone 15 32615 See map below
UTM, Zone 16 32616 See map below
UTM, Zone 17 32617 See map below
UTM, Zone 18 32618 See map below
UTM, Zone 19 32619 See map below

By Chrismurf at English Wikipedia, CC BY 3.0, https://commons.wikimedia.org/w/index.php?curid=40690482

Transforming to a different coordinate system

To transform our Lake Champlain shape file from NAD83 to WGS84, we will use the st_transform() function in the sf package. We will name the transformed shape file object lcsf.

lcsf<- lc %>% 
  st_transform("WGS84")

Run the lcsf object and plot it to confirm that the transformation was successful and that the geometry of the shape looks correct.

# Print lcsf shape file object
lcsf
## Simple feature collection with 1 feature and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.48796 ymin: 43.53135 xmax: -73.07608 ymax: 45.08523
## Geodetic CRS:  WGS 84
##   OBJECTID COMID    GNIS_NAME AREASQKM      REACHCODE    FTYPE GDB_GEOMAT
## 1    29848  <NA> Narrows, The 1146.091 02010008009170 LakePond       <NA>
##   SHAPESTAre SHAPESTLen                       geometry
## 1 1141293318    1102278 POLYGON ((-73.14224 45.0851...
# Plot lcsf
ggplot(lcsf) +
  geom_sf()

Now that our shape file is in the right coordinate system, we can add it to our map of New York and Vermont by adding an additional geom_sf() layer.

ggplot(ny_vt) +
  # add state outlines
  geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
  # add lake champlain shape
  geom_sf(data = lcsf, color = "#1282A2", fill = "#1282A2") +
  # remove background
  theme_void()

Setting a coordinate system

So transforming the Lake Champlain shape file to the correct coordinate system was fairly straightforward. Let’s take a look at a trickier case. We’ll import the shape file that represents the extent of the 2018 forest fire at the Altona Flat Rock pine barrens.

First, we’ll import the shapefile from the mapfiles folder using the st_read() function.

# Import shape file data from mapfiles folder
fr <- st_read("mapfiles", "fire_perimeter")
## Reading layer `fire_perimeter' from data source 
##   `C:\Users\malld001\Dropbox\PlattsburghCourses\ENV422_EnvDataAnalysis_Spring2024\A10Mapping\mapfiles' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -8198195 ymin: 5599219 xmax: -8195719 ymax: 5602464
## CRS:           NA

Note that no coordinate system is listed under CRS. Moreover, the coordinates for the bounding box look very different from the latitudes and longitudes we’ve seen so far. That means we’re dealing with a very different coordinate system, and we don’t know what it is! We will need to use the st_set_crs() function to specify the current coordinate system before we can transform the data to the “WGS84” system that we are using. The syntax for this function is
dataframe %>% st_set_crs(EPSG)

After much trial and error (set -> transform -> plot -> work?), I determined that the coordinates were in Mercator (EPSG = 3857). We will use st_set_crs() to set the object to its current projection and then st_transform() transform it to WGS84. We’ll save the transformed shape file as frsf.

frsf <- fr %>% 
  st_set_crs(3857) %>% 
  st_transform("WGS84")
frsf
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.64564 ymin: 44.85816 xmax: -73.62339 ymax: 44.87882
## Geodetic CRS:  WGS 84
##   Id                       geometry
## 1  0 POLYGON ((-73.64564 44.8639...

Note that the new frsf object has a coordinate system assigned, and that the Bounding box now looks like it contains latitude and longitudes in a sensible region of the planet. Try plotting the shape file on its own and on the ny_vt map to see how it looks at the different scales!

Combining map layers

You can also combine layers available in the maps package by adding geom_sf() layers just as we did for the shape file.

In this section, we will create a map of New York with the boundaries of counties included. First we will use map() to import the county polygons and st_as_sf() to transform the data to a spatial data frame, just as we did for states.

# Get counties from maps package
county <- maps::map("county", fill = T, plot = F)

# Create sf object
countysf <- st_as_sf(county)

Take a look at the countysf object. Notice that the ID contains a single column with both the county and state names. As we learned in the data management section, this format is not ideal. However, we can use functions in the dplyr package to easily separate this column into two columns, one for state and one for county.

# Overwrite data frame with new one that contains separate state and county columns
countysf <- countysf %>% 
  # create a new column, ID2, that is identical to ID
  mutate(ID2 = ID) %>% 
  # separate the new column into two columns, using "," as the delimiter
  separate(ID2, into = c("state", "county"), sep = ",")

When you view the countysf data frame now, you will see the original column as well as your new separated columns containing the state and county names.

Now let’s use the tidyverse to select only counties in New York and create a map.

countysf %>% 
  # select only counties in new york
  filter(state %in% c("new york")) %>% 
  # create plot and add geom_sf layer
  ggplot() +
  geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
  # eliminate background features
  theme_void()

Once we have this map, we can add as many layers as we like. As an example, let’s add a layer with the outline of New York.

# Create a data frame with just the outline of New York from the statesf data frame
ny <- statesf %>% 
  filter(ID %in% "new york")

# Create a data frame that contains only the county boundaries in NY from the countysf data frame
county_ny <- countysf %>% 
  filter(state %in% c("new york"))

# Copy map code from above
ggplot(county_ny) +
  geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
  theme_void() +
  # add layer for outline of New York state
  geom_sf(data = ny, color = "#283044", fill = NA)

Challenge 2

Based on what you have learned, create a map for the state of New York that includes county boundaries and a polygon showing the outline of Lake Champlain.
Add a layer that places a labelled point over SUNY Plattsburgh.
When you are satisfied with your map, use ggsave("A10_C2.jpg") to export your map image.
Upload your plot to Challenge 2 on Brightspace.

Answer 
# Create an object that contains only the county boundaries in NY
county_ny <- countysf %>% 
  filter(state %in% c("new york"))
# Create map
ggplot(county_ny) +
  # add county boundaries
  geom_sf(fill = "#74D3AE", color = "#F5DBCB") +
  theme_void() +
  # add layer for outline of New York state
  geom_sf(data = ny, color = "#283044", fill = NA) +
  # add layer for outline of Lake Champlain
  geom_sf(data = lcsf, color = "#1282A2", fill = "#1282A2") +
  # add point at SUNY Plattsburgh coordinates
  geom_point(aes(x = -73.464279, y = 44.692614), size = 3, color = "#283044") +
  # add text at SUNY Plattsburgh coordinates and adjust position
  geom_text(aes(x = -73.464279, y = 44.692614, label = "SUNY Plattsburgh"),
            hjust = 1, nudge_x = -0.2, color = "#283044") +
  # remove background elements
  theme_void()

ggsave("figures/A10_C2.jpg")
## Saving 7 x 5 in image

Adding Data

Now that we can create a nice map of the counties, ideally we would also like to add data. To do this, we will first import a data frame that I downloaded from Wikipedia and cleaned in Excel. We will import this file from the mapfiles folder that you downloaded from Brightspace. When you view the data frame, note that I added a column for “county” that has the name of the counties exactly as they appear in the county_ny data frame.

#Import wikipedia data
ny_wiki <- read_csv("mapfiles/nycounty_info.csv")
## Rows: 62 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): county, CountyName, CountySeat, formed_from, named_for
## dbl (2): FIPS_Code, est
## num (2): pop_2010, Area_km2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#View first few lines of the data frame
head(ny_wiki)
## # A tibble: 6 × 9
##   county      CountyName  FIPS_C…¹ Count…²   est forme…³ named…⁴ pop_2…⁵ Area_…⁶
##   <chr>       <chr>          <dbl> <chr>   <dbl> <chr>   <chr>     <dbl>   <dbl>
## 1 albany      Albany             1 Albany   1683 "One o… "James…  304204    1380
## 2 allegany    Allegany           3 Belmont  1806 "Genes… "A var…   48946    2678
## 3 bronx       Bronx              5 none     1914 "New Y… "Jonas… 1385108     149
## 4 broome      Broome             7 Bingha…  1806 "Tioga… "John …  200600    1852
## 5 cattaraugus Cattaraugus        9 Little…  1808 "Genes… "A wor…   80317    3393
## 6 cayuga      Cayuga            11 Auburn   1799 "Onond… "The\x…   80026    2238
## # … with abbreviated variable names ¹​FIPS_Code, ²​CountySeat, ³​formed_from,
## #   ⁴​named_for, ⁵​pop_2010, ⁶​Area_km2

Now we will join our ny_wiki information data frame to our county_ny data frame using a left_join() from the dplyr package. (Consult your dplyr cheat sheet for more information on available join functions. You may find these useful as you continue to work with your project data!) We will also use the mutate() function to calculate a people_per_km2 variable to express population density.

# Create joined data frame
ny_info <- left_join(county_ny, ny_wiki, by = "county") %>% 
  # Calculate population density
  mutate(people_per_km = pop_2010 / Area_km2)
# Use str() or head() to inspect your new data frame

Plot population density by county

Using our new dataframe ny_info, we will plot the map of New York by referencing our basemap and adding a layer that fills county boundaries based on our new people_per_km variable.

#Create new map object with layer for pop density
ggplot(ny_info) +
  # add sf geometry with fill corresponding to population density
  geom_sf(aes(fill = people_per_km), color = "#F5DBCB") +
  theme_void()

Hmm…all the counties look pretty much the same. The problem is that the population of New York City is so “yuge” that it makes it hard to see differences among the other counties. We can fix this problem by rescaling the gradient. This is similar to performing a log10 transformation on the y-axis to better observe differences among groups.

#Create new map object with layer for pop density
ggplot(ny_info) +
  # add sf geometry with fill corresponding to population density
  geom_sf(aes(fill = people_per_km), color = "#F5DBCB") +
  theme_void() +
  # add layer to density map rescaling the gradient by log10
  scale_fill_gradient(trans = "log10")

Challenge 3

Based on what you have learned, create a map for the state of New York that color codes counties based on the year they were established.
Consult Google or your ggplot cheat sheet to learn how to change the color gradient for your map.
Consult Google or your ggplot cheat sheet to learn how to change the title of your legend.
When you are satisfied with your map, use ggsave("A10_C3.jpg") to export your map image.
Upload your plot to Challenge 3 on Brightspace.

Answer 
#Create new map object with layer for time established
ggplot(ny_info) +
  # add sf geometry with fill corresponding to year of establishment
  geom_sf(aes(fill = est), color = "white") +
  # change color of gradient and name of legend
  scale_fill_gradient(low = "#F5DBCB", high = "#8B5666", name = "Year Established") +
  theme_void()

#Export image to figures folder
ggsave("figures/A10_C3.jpg")
## Saving 7 x 5 in image

Zooming in on a map

Let’s try zooming in on our population density map to just look at the area around Clinton County. When I view this sf object, it will give me coordinates for the “Bounding box” which I can use to set limits on the scale of my graph. I will switch back to theme_minimal() so that we can see the latitude and longitude.

# Create data frame with rows of only Clinton county
clinton <- ny_info[ny_info$county == "clinton",]
# Then I will use the information to zoom into my population density map
# Focus on coordinates of the Bounding box
clinton
## Simple feature collection with 1 feature and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.03188 ymin: 44.43861 xmax: -73.34433 ymax: 45.01157
## Geodetic CRS:  +proj=longlat +ellps=clrk66 +no_defs +type=crs
##                  ID    state  county CountyName FIPS_Code  CountySeat  est
## 10 new york,clinton new york clinton    Clinton        19 Plattsburgh 1788
##          formed_from
## 10 Washington County
##                                                                                                                         named_for
## 10 George Clinton\xa0(1739\x961812), fourth\xa0Vice President of the United States\xa0and first and third\xa0Governor of New York
##    pop_2010 Area_km2                           geom people_per_km
## 10    82128     2896 MULTIPOLYGON (((-74.03188 4...      28.35912

I’m going to use the boundary of this sf object as a guide to zoom into my population density map using the function coord_sf().

ggplot(ny_info) +
  geom_sf(aes(fill = people_per_km), color = "#F5DBCB") +
  theme_minimal() +
  scale_fill_gradient(trans = "log10") +
  coord_sf(xlim = c(-74.2, -73.3), ylim = c(44.3, 45.1))

Note that I definitely had to tweak the xlim and ylim a bit to get the map to look the way I want. Being able to see the axes makes this a bit easier to do. You can always go back and add theme_void() to remove the axes once you have adjusted the limits.

Challenge 4

Based on what you have learned, create a map of population density that zooms into the area around Manhattan and Queens. (Hint: Manhattan is in “New York” County, and Queens is in “Queens” County.)
Use theme_minimal() to view and adjust your xlim and ylim as needed.
When you are satisfied with your map, use ggsave("A10_C4.jpg") to export your map image.
Upload your plot to Challenge 4 on Brightspace.

Answer 
#Find reasonable max and min starting points
nyc <- ny_info[ny_info$county == "new york"| ny_info$county == "queens",]
nyc
## Simple feature collection with 2 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.02615 ymin: 40.54823 xmax: -73.71102 ymax: 40.84044
## Geodetic CRS:  +proj=longlat +ellps=clrk66 +no_defs +type=crs
##                   ID    state   county CountyName FIPS_Code CountySeat  est
## 31 new york,new york new york new york   New York        61       none 1683
## 41   new york,queens new york   queens     Queens        81       none 1683
##                                                      formed_from
## 31 One of 12 original counties created in the\xa0New York colony
## 41 One of 12 original counties created in the\xa0New York colony
##                                                                                                                                                        named_for
## 31 King\xa0James II of England\xa0(1633\x961701), who was Duke of York and Albany before he ascended the throne of England, Duke of York being his English title
## 41                                                           Catherine of Braganza\xa0(1638\x961705), Queen of England and wife of King\xa0Charles II of England
##    pop_2010 Area_km2                           geom people_per_km
## 31  1585873       87 MULTIPOLYGON (((-73.92874 4...     18228.425
## 41  2230722      462 MULTIPOLYGON (((-73.76259 4...      4828.403
#Plot population density
ggplot(ny_info) +
  geom_sf(aes(fill = people_per_km), color = "#F5DBCB") +
  theme_minimal() +
  scale_fill_gradient(trans = "log10") +
  coord_sf(xlim = c(-74.05 , -73.5 ), ylim = c(40.5 , 40.9))

ggsave("figures/A10_C4.jpg")
## Saving 7 x 5 in image

You will notice that this map does not look great. The built-in maps often do not do well around coastlines when you zoom in. For coastline maps, you would be better served by data resources such as NOAA’s GSHHS: A Global Self-consistent, Hierarchical, High-resolution Geography Database.
https://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html

Importing and plotting raster data

So far, all of our mapping has involved vector graphics, or data based on lines that connect points in x-y space. Mapping data may also contain raster data, or data stored as pixels at various resolutions. These are the same types of data stored in photographs, and are commonly used in remote sensing and in various types of aerial imagery. Here we will use the terrainr and raster packages to plot aerial imagery of the Lake Champlain Basin and the Altona Flat Rock burn site. This code should work for any georeferenced aerial imagery stored as a .tif. Note that I also used the terrainr package to download this aerial imagery from the USGS database. For more information on how to import imagery, see the complete introduction to terrainr in Mahoney (2022).

First, we will use the raster::stack() function to import the data for the Lake Champlain Basin. Note that this is a very large geographic area, so the data were imported at a resolution of 100 m. This means that each pixel on the map represents a 100 m x 100 m area of land. Then, we will convert the pixel data to a data frame object using the as.data.frame() function, and define the names of the columns using the setNames() function.

# Import raster data from .tif file
lc_raster <- raster::stack("mapfiles/lakechamp_ortho_USGS.tif")
# Convert raster data 
lc_raster_df <- as.data.frame(lc_raster, xy = TRUE)
lc_raster_df <- setNames(lc_raster_df, c("x", "y", "red", "green", "blue"))

Take a closer look at the lc_raster_df data frame. It contains x and y data, which represent longitudes and latitudes of the center of each pixel, as well as rgb color values that specify the color of the pixel. Note that these are stored as numeric values, so you can analyze these or perform calculations as you would any other numeric variable. Such analysis is beyond the scope of this class but forms the basis of most remote sensing. Let’s plot the data as a map layer using geom_spatial_rgb() function in the terrainr package. Note that we also add the coord_sf("WGS84") function from the sf package to specify the projection we want to use. (This plot contains a lot of information, so it will take a few seconds to generate!)

ggplot() + 
  # Add raster layer
  geom_spatial_rgb(data = lc_raster_df,
                   # Required aesthetics r/g/b specify color bands:
                   aes(x = x, y = y, r = red, g = green, b = blue)) + 
  # Specify coordinate system
  coord_sf(crs = "WGS84")

So that’s a great start! Notice that the area in Canada projects as missing data, which is plotted as a rather unflattering shade of pink. Let’s adjust the zoom on the map to exclude Canada by adding a ylim = argument to our coord_sf() function. While we’re at it, we’ll add the Lake Champlain shape file back to the map as an outline and remove the x and y axis labels.

ggplot() + 
  # Add raster layer
  geom_spatial_rgb(data = lc_raster_df,
                   # Required aesthetics r/g/b specify color bands:
                   aes(x = x, y = y, r = red, g = green, b = blue)) + 
  # Add lcsf shape file
  geom_sf(data = lcsf, color = "#1282A2", fill = NA) +
  # Change longitude limits
  coord_sf(crs = 4326, ylim = c(43.3, 44.95)) +
  # Remove background shading
  theme_minimal() +
  # Remove x and y labels
  labs(x = "",
       y = "")

Adding scale bars and north arrows

Scale bars are essential if one of the goals of your map is to allow a user to determine distances between locations. Similarly, a north arrow should always be included to indicate the orientation of the map. The ggspatial package allows you to add a scale bar and north arrow to any ggplot mapping object quickly, and it has several tools that allow you to easily customize the scale and appearance of these objects. Let’s add a scale bar to our Lake Champlain app using the annotation_scale() function, and a north arrow using the annotation_north_arrow() function. Note that both of these functions will default to placing the element in the bottom left corner of the map, so we will add the location = "tl" argument to our north arrow function to specify that we would like it to go in the top left of the map.

ggplot() + 
  # Add raster layer
  geom_spatial_rgb(data = lc_raster_df,
                   # Required aesthetics r/g/b specify color bands:
                   aes(x = x, y = y, r = red, g = green, b = blue)) + 
  # Add lcsf shape file
  geom_sf(data = lcsf, color = "#1282A2", fill = NA) +
  # Change longitude limits
  coord_sf(crs = 4326, ylim = c(43.3, 44.95)) +
  # Remove background shading
  theme_minimal() +
  # Remove x and y labels
  labs(x = "",
       y = "") +
  # Add a scale bar using annotation_scale()
  annotation_scale() +
  # Add a north arrow using annotation_north_arrow(location = "tl")
  annotation_north_arrow(location = "tl")

Load the help file for the annotation functions to look at the customization options. (hint: run ?annotation_north_arrow in the console)

Challenge 5

Based on what you have learned, plot the aerial imagery for the Altona Flat Rock stored in the “flatrock_ortho_USGS.tif” file in the mapfiles folder.
Add the “frsf” shape file to the map to outline the perimeter of the 2018 burn, and change the color of the outline so that it’s easy to see.
Add a scale bar and north arrow to your map.
Note: You may need to adjust limits to y and x axis in the coord_sf() function for map to look right.
When you are satisfied with your map, use ggsave("A10_C5.jpg") to export your map image.
Upload your plot to Challenge 5 on Brightspace.

Answer 
# Import raster data from .tif file
fr_raster <- raster::stack("mapfiles/flatrock_ortho_USGS.tif")
# Convert raster data 
fr_raster_df <- as.data.frame(fr_raster, xy = TRUE)
fr_raster_df <- setNames(fr_raster_df, c("x", "y", "red", "green", "blue"))

ggplot() + 
  # Add raster layer
  geom_spatial_rgb(data = fr_raster_df,
                   # Required aesthetics r/g/b specify color bands:
                   aes(x = x, y = y, r = red, g = green, b = blue)) + 
  # Add frsf shape file
  geom_sf(data = frsf, color = "yellow", fill = NA) +
  # Change longitude limits
  coord_sf(crs = 4326, ylim = c(44.857, 44.88), xlim = c(-73.65, -73.62)) +
  # Remove background shading
  theme_minimal() +
  # Remove x and y labels
  labs(x = "",
       y = "") +
  # Add a scale bar using annotation_scale()
  annotation_scale() +
  # Add a north arrow using annotation_north_arrow(location = "tr")
  annotation_north_arrow(location = "tr")

ggsave("figures/A10_C5.jpg")
## Saving 7 x 5 in image

NEXT WEEK: Group Assignment

Based on what you have learned, create a map for your data project.
Consider adding: points and labels for plot locations, color-coding points for plot locations based on treatments, changing the size/color of points to correspond to variables in your dataset.
ENV 422: Please feel free to exchange ideas and code with your other group members.
When you are satisfied with your map, use ggsave() to export your map image and upload it to the Assignment 11 dropbox on Brightspace.

References

Anderson, Eric C., Kristen C. Ruegg, Tina Cheng, et al. 2017. Chapter 8. Making simple maps with R. Case studies in reproducible research: a spring seminar at UCSC. Online. https://eriqande.github.io/rep-res-eeb-2017/map-making-in-R.html

Dunnington, Dewey. 2021. ggspatial: Spatial data framework for ggplot2. Online. https://cran.r-project.org/web/packages/ggspatial/ggspatial.pdf

Frazier, Melanie. 2020. Overview of Coordinate Reference Systems (CRS) in R. Online. https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf

Lovelace, Robin, Jakub Nowosad, and Jannes Muenchow. 2021. Geocomputation with R. Open Source Textbook. CRC Press. https://bookdown.org/robinlovelace/geocompr/

Mahoney, Michael. 2022. A gentle introduction to terrainr. Online. https://cran.r-project.org/web/packages/terrainr/vignettes/overview.html